rm(list=ls())
setwd("Main Text-Figures/Figure6/")

#############################################################################################

# simulations for dimensionality result across levels of prevalence 

d.prev <- c(NA,NA,NA,NA)
yd <- c(NA,NA,NA,NA)
xd <- c(NA,NA,NA,NA)

  load("dimensionality/dimension-data-k35-v500-data.RData")
    d.prev[1] <- mean(sim.info.all$percent.rcv)
    yd[1] <- as.numeric(table(eigen.rcv[2,]<=1)[2])/(dim(eigen.rcv)[2])*100
    xd[1] <- mean(sim.info.all$rcv)
  
  load("dimensionality/dimension-data-k35-v750-data.RData")
    d.prev[2] <- mean(sim.info.all$percent.rcv)
    yd[2] <- as.numeric(table(eigen.rcv[2,]<=1)[2])/(dim(eigen.rcv)[2])*100
    xd[2] <- mean(sim.info.all$rcv)

  load("dimensionality/dimension-data-k35-v1000-data.RData")
    d.prev[3] <- mean(sim.info.all$percent.rcv)
    yd[3] <- as.numeric(table(eigen.rcv[2,]<=1)[2])/(dim(eigen.rcv)[2])*100
    xd[3] <- mean(sim.info.all$rcv)

  load("dimensionality/dimension-data-k35-v1500-data.RData")
    d.prev[4] <- mean(sim.info.all$percent.rcv)
    yd[4] <- as.numeric(table(eigen.rcv[2,]<=1)[2])/(dim(eigen.rcv)[2])*100
    xd[4] <- mean(sim.info.all$rcv)

#############################################################################################
    
# simulations for cohesion result across levels of prevalence 

co.prev <- c(NA,NA,NA,NA,NA)
yco <- c(NA,NA,NA,NA,NA)
xco <- c(NA,NA,NA,NA,NA)
  
  load("cohesion/cohesion-density-v500-data.RData")
    co.prev[1] <- mean(totals.dat$prevalence)*100
    yco[1] <- mean(totals.dat$codiff.over.10)
    xco[1] <- mean(totals.dat$rcv.total)

  load("cohesion/cohesion-density-v750-data.RData")
    co.prev[2] <- mean(totals.dat$prevalence)*100
    yco[2] <- mean(totals.dat$codiff.over.10)
    xco[2] <- mean(totals.dat$rcv.total)

  load("cohesion/cohesion-density-v1000-data.RData")
    co.prev[3] <- mean(totals.dat$prevalence)*100
    yco[3] <- mean(totals.dat$codiff.over.10)
    xco[3] <- mean(totals.dat$rcv.total)

  load("cohesion/cohesion-density-v1500-data.RData")
    co.prev[4] <- mean(totals.dat$prevalence)*100
    yco[4] <- mean(totals.dat$codiff.over.10)
    xco[4] <- mean(totals.dat$rcv.total)

  load("cohesion/cohesion-density-v2000-data.RData")
    co.prev[5] <- mean(totals.dat$prevalence)*100
    yco[5] <- mean(totals.dat$codiff.over.10)
    xco[5] <- mean(totals.dat$rcv.total)

#############################################################################################

# simulations for ideal point estimation result across levels of prevalence 
    
idp.prev <- c(NA,NA,NA,NA,NA)
yidp <- c(NA,NA,NA,NA,NA)
xidp <- c(NA,NA,NA,NA,NA)
    
  load("ideal-points/prev-idp-v500-data.RData")
    idp.prev[1] <- mean(prevalence)
    yidp[1] <- mean(total.median.rcv-total.median.full)
    xidp[1] <- mean(totalrcv)
    
  load("ideal-points/prev-idp-v750-data.RData")
    idp.prev[3] <- mean(prevalence)
    yidp[3] <- mean(total.median.rcv-total.median.full)
    xidp[3] <- mean(totalrcv)
    
  load("ideal-points/prev-idp-v1000-data.RData")
    idp.prev[2] <- mean(prevalence)
    yidp[2] <- mean(total.median.rcv-total.median.full)
    xidp[2] <- mean(totalrcv)

  load("ideal-points/prev-idp-v1500-data.RData")
    idp.prev[4] <- mean(prevalence)
    yidp[4] <- mean(total.median.rcv-total.median.full)
    xidp[4] <- mean(totalrcv)

  load("ideal-points/prev-idp-v2000-data.RData")
    idp.prev[5] <- mean(prevalence)
    yidp[5] <- mean(total.median.rcv-total.median.full)
    xidp[5] <- mean(totalrcv)
    
#############################################################################################

#Figure 6
    
  pdf("Fig6_.pdf",width=6,height=4)
    par(mar=c(4,6.5,1,3.5))
    plot(x=xd,y=yd,type='l',lty=2,lwd=2,ylim=c(0,52),xlim=c(40,210),
         ylab="",xlab="",xaxt='n',yaxt='n')
    points(x=xd,y=yd,pch=22,bg="white")
    mtext(side=2,line=3.5,"Percentage of Legislatures\n with RCV Inference Problem",cex=1.25)
    mtext(side=1,line=3,"Number of Roll Call Votes",cex=1.25)
    axis(side=1,at=c(50,75,100,125,150,175,200),labels=c(50,"",100,"",150,"",200))
    axis(side=2,at=seq(0,100,25),labels=c("0%","25%","50%","75%","100%"),las=1)
    lines(x=xco,y=yco,col="grey50",lwd=4)
    axis(side=4,at=c(0,10,25,50),labels=c("0%","10%","25%","50%"),las=1)
    points(x=xco,y=yco,pch=23,bg="grey60",col="grey60")
    lines(x=xidp,y=yidp,lwd=2)
    points(x=xidp,y=yidp,pch=24,bg="black")
    
    
    w <- 1.5
    for(i in 1:length(xco)){
      polygon(x=c(xco[i]-w,xco[i]-w,xco[i]+w,xco[i]+w),
              y=c(0,co.prev[i],co.prev[i],0),
              col="grey60",border=NA)
    }
    for(i in 1:length(xd)){
      polygon(x=c(xd[i]-w,xd[i]-w,xd[i]+w,xd[i]+w),
              y=c(0,d.prev[i],d.prev[i],0),
              col="white")
    }
    for(i in 1:length(xidp)){
      polygon(x=c(xidp[i]-w,xidp[i]-w,xidp[i]+w,xidp[i]+w),
              y=c(0,idp.prev[i],idp.prev[i],0),
              col="black",border=NA)
    }
    abline(h=0)
  dev.off()